######################################
## Acharya, Blackwell, Sen (2018)
## Replication of Sen (2017)


######################################
## Necessary packages
######################################

library(stargazer)
library(xtable)
library(reshape)
library(foreign)
library(lmtest)
library(sandwich)
library(ggplot2)


# function that does clustered SEs

robust.se <- function(fm, clvar){
  vcovCL <- vcovClustered(fm, clvar)
  coeftest(fm, vcovCL)
}

vcovClustered <- function(fm, clvar){
    # R-codes (www.r-project.org) for computing
    # clustered-standard errors. Mahmood Arai, Jan 26, 2008.
    # The arguments of the function are:
    # fitted model, cluster1 and cluster2
    # You need to install libraries `sandwich' and `lmtest'
  library(sandwich);library(lmtest);
  ##env <- environment(fm$terms)
  x <- eval(fm$call$data, envir = parent.frame())
  ##browser()
  if (class(fm) == "polr") {
    require(MASS)
    cluster <- x[rownames(predict(fm, type = "probs")), clvar]
  } else {
    cluster <- x[names(predict(fm)), clvar]
  }
  if (is.factor(cluster)) {
        cluster <- factor(cluster)
  }
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- dim(vcov(fm))[1]
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  return(vcovCL)
}


######################################
## Loading the data
######################################

conjoint <- read.csv("Judicial_Nominations_FINAL.csv",stringsAsFactors = FALSE)
conjoint <- conjoint[-1,]
dim(conjoint)

conjoint <- subset(conjoint, Q2.1 != "")


######################################
## Creating Some Variables for Analysis
######################################


West <- c("CA","NV","AZ","NM","OR","WA","HI","UT","AK","ID","WY","CO","MT")
North.East <- c("ME","NH","VT","CT","MA","RI","NY","NJ","PA","DE","MD","DC")
South <- c("VA","WV","TN","KY","AL","MS","GA","FL","SC","NC","OK","TX","LA")
Midwest <- c("KS","MO","OH","MI","IN","IL","NE","WI","MN","IA","AK","SD","ND")
length(West) + length(North.East) + length(South) + length(Midwest)

## creating a "Region" variable

conjoint$region <- NA
conjoint$region[which(conjoint$Q2.2 %in% West)] <- "West"
conjoint$region[which(conjoint$Q2.2 %in% South)] <- "South"
conjoint$region[which(conjoint$Q2.2 %in% North.East)] <- "North East"
conjoint$region[which(conjoint$Q2.2 %in% Midwest)] <- "Midwest"
data.frame(conjoint$region, conjoint$Q2.2)

## creating a composite "Knowledge" variable

conjoint$knowledge <- NA
conjoint$knowledge <- as.numeric(conjoint$Q3.3 == "Appointed to the bench")+as.numeric(conjoint$Q3.4 == "Life term") + as.numeric(conjoint$Q3.5 == "U.S. Supreme Court") + as.numeric(conjoint$Q3.6 == "California") + as.numeric(conjoint$Q3.7 == "John Roberts") + as.numeric(conjoint$Q3.8 == "Elena Kagan")
table(conjoint$knowledge)

## creating a "Legitimacy" variable

conjoint$legitimacy <- NA
conjoint$legitimacy <- as.numeric(conjoint$Q4.2 == "Disagree"|conjoint$Q4.2 == "Strongly Disagree")+as.numeric(conjoint$Q4.3 == "Agree"|conjoint$Q4.3 == "Strongly Agree") + as.numeric(conjoint$Q4.4 == "Disagree"|conjoint$Q4.4 == "Strongly Disagree") + as.numeric(conjoint$Q4.5 == "Disagree"|conjoint$Q4.5 == "Strongly Disagree") + as.numeric(conjoint$Q4.6 == "Disagree"|conjoint$Q4.6 == "Strongly Disagree")
table(conjoint$legitimacy)

## creating a more straightforward party variable

conjoint$respondent.party <- NA
conjoint$respondent.party[which(
        conjoint$Q2.5 == "A Republican"|
        conjoint$Q2.5 == "A Strong Republican"
        #|conjoint$Q2.5 == "An Independent who leans Republican"
)] <- "Republican"


conjoint$respondent.party[which(
        conjoint$Q2.5 == "A Democrat"|
        conjoint$Q2.5 == "A Strong Democrat"
        #|conjoint$Q2.5 == "An Independent who leans Democrat"
)] <- "Democrat"
conjoint$respondent.party[which(conjoint$Q2.5 == "An Independent")] <- "Independent"


######################################
## Creating Some Subsets for Analysis
######################################

all <- conjoint
republicans <- subset(conjoint, Q2.5 == "A Republican"| Q2.5 == "A Strong Republican"| Q2.5 == "An Independent who leans Republican")
dim(republicans)
democrats <- subset(conjoint, Q2.5 == "A Democrat"| Q2.5 == "A Strong Democrat"| Q2.5 == "An Independent who leans Democrat")
dim(democrats)
independents <- subset(conjoint, Q2.5 == "An Independent")
dim(independents)

women <- subset(conjoint, Q2.3 == "Female")
men <- subset(conjoint, Q2.3 == "Male")

whites <- subset(conjoint, Q2.4 == "White, non-Hispanic")
blacks <- subset(conjoint, Q2.4 == "Black/African-American, non-Hispanic")
latinos <- subset(conjoint, Q2.4 == "Hispanic, Latino/a, or Chicano/a")
asians <- subset(conjoint, Q2.4 == "Asian, non-Hispanic")

#conjoint <- democrats


demographics <- conjoint[,c("Q1.2","Q2.1","Q2.2","Q2.3",
                "Q2.4","Q2.5","Q2.6","Q2.7")]

knowledge <- conjoint[,c("Q3.1","Q3.2","Q3.3",
                "Q3.4","Q3.5","Q3.6","Q3.7","Q3.8")]

legitimacy <- conjoint[,c("Q4.1","Q4.2","Q4.3",
                "Q4.4","Q4.5","Q4.6")]

id <- conjoint[,"pid"]


######################################
## Conjoint Cleaning
######################################


conjoint.cleaning <- function(data, x1, x2, x3, x4, x5, x6, x7, x8 = NULL, answer1, answer2, answer3, demographics, knowledge, legitimacy, id){
        conjoint <- data

        if(!is.null(x8)){
        conjoint1 <- data.frame(matrix(data = NA, ncol = 8,
                nrow = nrow(conjoint)))
        colnames(conjoint1) <- c("Previous work experience","Gender",
                "Age","Race or ethnicity","Previous judicial clerkship experience",
                "Legal education","Religion","Party leaning")
        answer.cats <- cbind(conjoint[,answer1],
                conjoint[,answer2], conjoint[,answer3])
        colnames(answer.cats) <- c("Support","Qualified","Trust")
        }

        else{
        conjoint1 <- data.frame(matrix(data = NA, ncol = 7,
        nrow = nrow(conjoint)))
        colnames(conjoint1) <- c("Previous work experience","Gender",
                "Age","Race or ethnicity","Previous judicial clerkship experience",
                "Legal education","Religion")
        answer.cats <- cbind(conjoint[,answer1],
                conjoint[,answer2], conjoint[,answer3])
        colnames(answer.cats) <- c("Support","Qualified","Trust")
                }

        for(i in 1:ncol(conjoint1)){
                foo <- grep(names(conjoint1)[i], conjoint[,x1])
                bar <- conjoint[foo,which(colnames(conjoint) == x1)+1]
                conjoint1[foo,i] <- as.character(bar)
        }

        for(i in 1:ncol(conjoint1)){
                foo <- grep(names(conjoint1)[i], conjoint[,x2])
                bar <- conjoint[foo,which(colnames(conjoint) == x2)+1]
                conjoint1[foo,i] <- as.character(bar)
        }

        for(i in 1:ncol(conjoint1)){
                foo <- grep(names(conjoint1)[i], conjoint[,x3])
                bar <- conjoint[foo,which(colnames(conjoint) == x3)+1]
                conjoint1[foo,i] <- as.character(bar)
        }

        for(i in 1:ncol(conjoint1)){
                foo <- grep(names(conjoint1)[i], conjoint[,x4])
                bar <- conjoint[foo,which(colnames(conjoint) == x4)+1]
                conjoint1[foo,i] <- as.character(bar)
        }

        for(i in 1:ncol(conjoint1)){
                foo <- grep(names(conjoint1)[i], conjoint[,x5])
                bar <- conjoint[foo,which(colnames(conjoint) == x5)+1]
                conjoint1[foo,i] <- as.character(bar)
        }

        for(i in 1:ncol(conjoint1)){
                foo <- grep(names(conjoint1)[i], conjoint[,x6])
                bar <- conjoint[foo,which(colnames(conjoint) == x6)+1]
                conjoint1[foo,i] <- as.character(bar)
        }

        for(i in 1:ncol(conjoint1)){
                foo <- grep(names(conjoint1)[i], conjoint[,x7])
                bar <- conjoint[foo,which(colnames(conjoint) == x7)+1]
                conjoint1[foo,i] <- as.character(bar)
        }


        if(!is.null(x8)){
        for(i in 1:ncol(conjoint1)){
                foo <- grep(names(conjoint1)[i], conjoint[,x8])
                bar <- conjoint[foo,which(colnames(conjoint) == x8)+1]
                conjoint1[foo,i] <- as.character(bar)
        }

        conjoint1 <- cbind(conjoint1, answer.cats)
        colnames(conjoint1) <- c("Work","Gender","Age","Race","Clerk",
                "Educ","Religion","Party","Support","Qualification","Trust")
        }

        else{
        conjoint1 <- cbind(conjoint1, answer.cats)
        colnames(conjoint1) <- c("Work","Gender","Age","Race","Clerk",
                "Educ","Religion","Support","Qualification","Trust")
                }

        conjoint1$Trust <- factor(conjoint1$Trust,
                levels = c("Strongly mistrust","Mistrust",
                "Neither trust or mistrust",
                "Trust","Strongly trust"))
        conjoint1$Qualification <- factor(conjoint1$Qualification,
                levels = c("Highly unqualified","Unqualified",
                "Neither qualified nor unqualified","Qualified","Highly qualified"))
        conjoint1$Support <- factor(conjoint1$Support,
                levels = c("Strongly oppose","Oppose","Neither support nor oppose",
                "Support","Strongly support"))
        conjoint1$Race <- factor(conjoint1$Race,
                levels = c("White","Black","Hispanic or Latino/a","Asian American"))
        conjoint1$Work <- factor(conjoint1$Work,
                levels = c("Lawyer in private practice","Law professor",
                "Lower-level federal judge", "Non-profit lawyer",
                "Public defender","Prosecutor","State judge","Elected politician"))
        conjoint1$Religion <- factor(conjoint1$Religion,
                levels = c("Mainline Protestant","Catholic",
                "Jewish", "Evangelical Protestant",
                "Mormon"))
        conjoint1$Educ <- factor(conjoint1$Educ,
                levels = c("Law School Ranked in Top 14","Law School Ranked #15-25","Law School Ranked #25-50","Law School Ranked #50-100","Law School Not Ranked in Top 100"))

        return(cbind(id, conjoint1, demographics, knowledge, legitimacy))
}

conjoint.P1 <- conjoint.cleaning(x1 = "P.1.1", x2 = "P.1.2", x3 = "P.1.3",
        x4 = "P.1.4", x5 = "P.1.5", x6 = "P.1.6", x7 =  "P.1.7", x8 = "P.1.8",
        answer1 = "Q16.2",answer2 = "Q16.3", answer3 = "Q16.4", data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)

head(conjoint.P1)

conjoint.P2 <- conjoint.cleaning(x1 = "P.2.1", x2 = "P.2.2", x3 = "P.2.3",
        x4 = "P.2.4", x5 = "P.2.5", x6 = "P.2.6", x7 = "P.2.7", x8 = "P.2.8",
        answer1 = "Q17.2",answer2 = "Q17.3", answer3 = "Q17.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.P3 <- conjoint.cleaning(x1 = "P.3.1", x2 = "P.3.2", x3 = "P.3.3",
        x4 = "P.3.4", x5 = "P.3.5", x6 = "P.3.6", x7 = "P.3.7", x8 = "P.3.8",
        answer1 = "Q18.2",answer2 = "Q18.3", answer3 = "Q18.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.P4 <- conjoint.cleaning(x1 = "P.4.1", x2 = "P.4.2", x3 = "P.4.3",
        x4 = "P.4.4", x5 = "P.4.5", x6 = "P.4.6", x7 = "P.4.7", x8 = "P.4.8",
        answer1 = "Q19.2",answer2 = "Q19.3", answer3 = "Q19.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.P5 <- conjoint.cleaning(x1 = "P.5.1", x2 = "P.5.2", x3 = "P.5.3",
        x4 = "P.5.4", x5 = "P.5.5", x6 = "P.5.6", x7 = "P.5.7", x8 = "P.5.8",
        answer1 = "Q20.2",answer2 = "Q20.3", answer3 = "Q20.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.P6 <- conjoint.cleaning(x1 = "P.6.1", x2 = "P.6.2", x3 = "P.6.3",
        x4 = "P.6.4", x5 = "P.6.5", x6 = "P.6.6", x7 = "P.6.7", x8 = "P.6.8",
        answer1 = "Q21.2",answer2 = "Q21.3", answer3 = "Q21.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.partisan.data <- rbind(conjoint.P1,conjoint.P2,conjoint.P3,conjoint.P4,conjoint.P5,conjoint.P6)

conjoint.NP1 <- conjoint.cleaning(x1 = "NP.1.1", x2 = "NP.1.2", x3 = "NP.1.3",
        x4 = "NP.1.4", x5 = "NP.1.5", x6 = "NP.1.6", x7 = "NP.1.7",
        answer1 = "Q10.2",answer2 = "Q10.3", answer3 = "Q10.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.NP2 <- conjoint.cleaning(x1 = "NP.2.1", x2 = "NP.2.2", x3 = "NP.2.3",
        x4 = "NP.2.4", x5 = "NP.2.5", x6 = "NP.2.6", x7 = "NP.2.7",
        answer1 = "Q11.2",answer2 = "Q11.3", answer3 = "Q11.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.NP3 <- conjoint.cleaning(x1 = "NP.3.1", x2 = "NP.3.2", x3 = "NP.3.3",
        x4 = "NP.3.4", x5 = "NP.3.5", x6 = "NP.3.6", x7 = "NP.3.7",
        answer1 = "Q12.2",answer2 = "Q12.3", answer3 = "Q12.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.NP4 <- conjoint.cleaning(x1 = "NP.4.1", x2 = "NP.4.2", x3 = "NP.4.3",
        x4 = "NP.4.4", x5 = "NP.4.5", x6 = "NP.4.6", x7 = "NP.4.7",
        answer1 = "Q13.2",answer2 = "Q13.3", answer3 = "Q13.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.NP5 <- conjoint.cleaning(x1 = "NP.5.1", x2 = "NP.5.2", x3 = "NP.5.3",
        x4 = "NP.5.4", x5 = "NP.5.5", x6 = "NP.5.6", x7 = "NP.5.7",
        answer1 = "Q14.2",answer2 = "Q14.3", answer3 = "Q14.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.NP6 <- conjoint.cleaning(x1 = "NP.6.1", x2 = "NP.6.2", x3 = "NP.6.3",
        x4 = "NP.6.4", x5 = "NP.6.5", x6 = "NP.6.6", x7 = "NP.6.7",
        answer1 = "Q15.2",answer2 = "Q15.3", answer3 = "Q15.4",data = conjoint,
        demographics = demographics, knowledge = knowledge, legitimacy = legitimacy, id = id)


conjoint.nonpartisan.data <- rbind(conjoint.NP1,conjoint.NP2,conjoint.NP3,conjoint.NP4,conjoint.NP5,conjoint.NP6)

conjoint.nonpartisan.data <- conjoint.nonpartisan.data[,1:11]
conjoint.partisan.data <- conjoint.partisan.data[,1:12]
conjoint.vars <- conjoint[,420:ncol(conjoint)]

######################################
## Re-merging the conjoint experiments with the other data
######################################

conjoint.partisan.data$id <- as.character(conjoint.partisan.data$id)
partisan.data <- merge(conjoint.partisan.data, conjoint.vars, by.x = "id", by.y = "pid")
dim(conjoint.partisan.data)
dim(partisan.data)

conjoint.nonpartisan.data$id <- as.character(conjoint.nonpartisan.data$id)
nonpartisan.data <- merge(conjoint.nonpartisan.data, conjoint.vars, by.x = "id", by.y = "pid")
dim(conjoint.nonpartisan.data)
dim(nonpartisan.data)


partisan.data$party.treatment <- 1
partisan.data$party <- partisan.data$Party
partisan.data$Party <- NULL
nonpartisan.data$party.treatment  <- 0
nonpartisan.data$party  <- NA

partisan.data$copartisan <- with(partisan.data, 1*((party %in% c("Strong Republican", "Leans Republican")) & respondent.party == "Republican") + 1*((party %in% c("Strong Democrat", "Leans Democrat")) & respondent.party == "Democrat"))
nonpartisan.data$copartisan <- NA

nonpartisan.data$Female <- as.numeric((nonpartisan.data$Gender == "Female"))
partisan.data$Female <- as.numeric((partisan.data$Gender == "Female"))

partisan.data$Age.f <- factor(partisan.data$Age)
partisan.data$Clerk.f <- factor(partisan.data$Clerk)
partisan.data$party.f <- factor(partisan.data$party, levels = c("Strong Democrat", "Leans Democrat", "Independent", "Leans Republican", "Strong Republican"))
partisan.data$copartisan.f <- factor(partisan.data$copartisan)
partisan.data$Trust2 <- 1*(as.numeric(partisan.data$Trust) > 3)
partisan.data$Support2 <- 1*(as.numeric(partisan.data$Support) > 3)
partisan.data$Qualification2 <- 1*(as.numeric(partisan.data$Qualification) > 3)



nonpartisan.data$Age.f <- factor(nonpartisan.data$Age)
nonpartisan.data$Clerk.f <- factor(nonpartisan.data$Clerk)
nonpartisan.data$party.f <- factor(nonpartisan.data$party, levels = c("Strong Democrat", "Leans Democrat", "Independent", "Leans Republican", "Strong Republican"))
nonpartisan.data$copartisan.f <- factor(nonpartisan.data$copartisan)
nonpartisan.data$Trust2 <- 1*(as.numeric(nonpartisan.data$Trust) > 3)
nonpartisan.data$Support2 <- 1*(as.numeric(nonpartisan.data$Support) > 3)
nonpartisan.data$Qualification2 <- 1*(as.numeric(nonpartisan.data$Qualification) > 3)
## basic results, all respondents

all.data <- rbind(nonpartisan.data, partisan.data)



covs <- c("Female", "Race", "Age.f", "Religion", "Clerk.f", "Educ", "Work")
cov.levs <- c("Gender:     ", "Male", "Female", NA,
                  "Race:       ", "White", "Black", "Hispanic", "Asian", NA,
                  "Age:        ", "30-40", "40-50", "50-60", "60-70", "70+", NA,
                  "Religion:   ", "Mainline Protestant", "Catholic", "Jewish", "Evangelical", "Mormon", NA,

                  "Served as Clerk:", "No", "Yes", NA,
                  "Law School Rank:", "Top 14", "15-24" ,"25-49","50-100", "100+", NA,
              "Previous Work:  ", "Private Practice", "Law Professor",
              "Federal Judge", "Non-profit", "Public Defender", "Prosecutor",
              "State Judge", "Elected Politician")

all.data <- all.data[!is.na(all.data$Support),]
all.data$Trust2 <- 1*(as.numeric(all.data$Trust) >= 3)


all.data$Qualification2 <- 1*(as.numeric(all.data$Qualification) >= 3)
        ## Qualified or Highly Qualified

all.data$Support2 <- 1*(as.numeric(all.data$Support) >= 3)

all.data$party.treatment <- 1-all.data$party.treatment

dems.data <- subset(all.data, respondent.party == "Democrat")
dems.co.data <- subset(dems.data, copartisan == 1 | party.treatment == 1)
dems.nc.data <- subset(dems.data, copartisan == 0 | party.treatment == 1)
amce.data <- subset(dems.data, party.treatment == 1)

dim(dems.data)

supp.fun <- function(x, data, int = TRUE) {
  if (int) {
    form <- substitute(Support2 ~ i * party.treatment, list(i = as.name(x)))
  } else {
    form <- substitute(Support2 ~ i, list(i = as.name(x)))
  }
  mod <- lm(form, data = data)
  vcovCL <- vcovClustered(mod, "id")
  coefs <- coef(mod)
  crse <- sqrt(diag(vcovCL))
  cvals95 <- qt(0.975, df = mod$df.residual)
  ci95 <- matrix(NA, nrow = 2, ncol = length(coefs))
  ci95[1,] <- coefs - cvals95 * crse
  ci95[2,] <- coefs + cvals95 * crse
  cvals90 <- qt(0.95, df = mod$df.residual)
  ci90 <- matrix(NA, nrow = 2, ncol = length(coefs))
  ci90[1,] <- coefs - cvals90 * crse
  ci90[2,] <- coefs + cvals90 * crse
  cvals80 <- qt(0.90, df = mod$df.residual)
  ci80 <- matrix(NA, nrow = 2, ncol = length(coefs))
  ci80[1,] <- coefs - cvals80 * crse
  ci80[2,] <- coefs + cvals80 * crse
  return(list(coef=coefs, crse=crse, ci95 = ci95, ci90 = ci90, ci80 = ci80))
}

pad.qoi <- function(x, i) c(NA, x[[i]], NA)
pad.ci <- function(x, i) cbind(NA, x[[i]], NA)
support.mods <- lapply(covs, supp.fun, data = amce.data, int = FALSE)
support.nc.mods <- lapply(covs, supp.fun, data = dems.nc.data)
support.co.mods <- lapply(covs, supp.fun, data = dems.co.data)
supp.acmes <- unlist(lapply(support.mods, pad.qoi, i = "coef"))
supp.nc.coefs <- unlist(lapply(support.nc.mods, pad.qoi, i = "coef"))
supp.co.coefs <- unlist(lapply(support.co.mods, pad.qoi, i = "coef"))
ints <- grepl("party.treatment", names(supp.nc.coefs))
supp.acmes <- supp.acmes
supp.nc.diffs <- supp.nc.coefs[ints | is.na(supp.nc.coefs)]
supp.co.diffs <- supp.co.coefs[ints | is.na(supp.co.coefs)]
supp.acme.ci95 <- sapply(support.mods, pad.ci, i = "ci95")
supp.acme.ci95 <- do.call(cbind, supp.acme.ci95)
supp.acme.ci95 <- supp.acme.ci95
supp.acme.ci90 <- sapply(support.mods, pad.ci, i = "ci90")
supp.acme.ci90 <- do.call(cbind, supp.acme.ci90)
supp.acme.ci90 <- supp.acme.ci90
supp.acme.ci80 <- sapply(support.mods, pad.ci, i = "ci80")
supp.acme.ci80 <- do.call(cbind, supp.acme.ci80)
supp.acme.ci80 <- supp.acme.ci80
supp.nc.diff.ci95 <- sapply(support.nc.mods, pad.ci, i = "ci95")
supp.nc.diff.ci95 <- do.call(cbind, supp.nc.diff.ci95)
supp.nc.diff.ci95 <- supp.nc.diff.ci95[, ints | is.na(supp.nc.diff.ci95[1,])]
supp.nc.diff.ci90 <- sapply(support.nc.mods, pad.ci, i = "ci90")
supp.nc.diff.ci90 <- do.call(cbind, supp.nc.diff.ci90)
supp.nc.diff.ci90 <- supp.nc.diff.ci90[, ints | is.na(supp.nc.diff.ci90[1,])]
supp.nc.diff.ci80 <- sapply(support.nc.mods, pad.ci, i = "ci80")
supp.nc.diff.ci80 <- do.call(cbind, supp.nc.diff.ci80)
supp.nc.diff.ci80 <- supp.nc.diff.ci80[, ints | is.na(supp.nc.diff.ci80[1,])]
supp.co.diff.ci95 <- sapply(support.co.mods, pad.ci, i = "ci95")
supp.co.diff.ci95 <- do.call(cbind, supp.co.diff.ci95)
supp.co.diff.ci95 <- supp.co.diff.ci95[, ints | is.na(supp.co.diff.ci95[1,])]
supp.co.diff.ci90 <- sapply(support.co.mods, pad.ci, i = "ci90")
supp.co.diff.ci90 <- do.call(cbind, supp.co.diff.ci90)
supp.co.diff.ci90 <- supp.co.diff.ci90[, ints | is.na(supp.co.diff.ci90[1,])]
supp.co.diff.ci80 <- sapply(support.co.mods, pad.ci, i = "ci80")
supp.co.diff.ci80 <- do.call(cbind, supp.co.diff.ci80)
supp.co.diff.ci80 <- supp.co.diff.ci80[, ints | is.na(supp.co.diff.ci80[1,])]
## Set baselines to zero
supp.acmes[grep("(Intercept)", names(supp.acmes))] <- 0
supp.acme.ci95[, grep("(Intercept)", names(supp.acmes))] <- 0
supp.acme.ci90[, grep("(Intercept)", names(supp.acmes))] <- 0
supp.acme.ci80[, grep("(Intercept)", names(supp.acmes))] <- 0
supp.nc.diffs[grep("^party.treatment", names(supp.nc.diffs))] <- 0
supp.nc.diff.ci95[, grep("^party.treatment", names(supp.nc.diffs))] <- 0
supp.nc.diff.ci90[, grep("^party.treatment", names(supp.nc.diffs))] <- 0
supp.nc.diff.ci80[, grep("^party.treatment", names(supp.nc.diffs))] <- 0
supp.co.diffs[grep("^party.treatment", names(supp.co.diffs))] <- 0
supp.co.diff.ci95[, grep("^party.treatment", names(supp.co.diffs))] <- 0
supp.co.diff.ci90[, grep("^party.treatment", names(supp.co.diffs))] <- 0
supp.co.diff.ci80[, grep("^party.treatment", names(supp.co.diffs))] <- 0
names(supp.acmes) <- cov.levs
names(supp.nc.diffs) <- cov.levs
names(supp.co.diffs) <- cov.levs

cairo_pdf(file="conjoint-plot.pdf", height = 6.5, width = 12, pointsize = 9)
par(mar = c(4.5, 10, 2.1, 0.1), mfrow = c(1,3))
plot(x = NULL, y = NULL, xlim = c(-0.2, 0.2), ylim = c(1,length(supp.acmes)),
     pch = 19, yaxt = "n", bty = "n", xlab = "Marginal Effects, Support",
     ylab = "", main = "Overall Marginal Effects")
abline(h = 1:length(supp.acmes), col = "grey90")
abline(v = 0, col = "grey60")
points(x = supp.acmes, y = length(supp.acmes):1, pch = 19)
segments(x0 = supp.acme.ci95[1,], x1 = supp.acme.ci95[2,],
         y0 = length(supp.acmes):1)
segments(x0 = supp.acme.ci90[1,], x1 = supp.acme.ci90[2,],
         y0 = length(supp.acmes):1, lwd = 2)
axis(side = 2, at = 1:length(supp.acmes),
     labels = rep("", times = length(supp.acmes)))
mtext(side = 2, at = length(supp.acmes):1, names(supp.acmes), las = 1,
      cex = 0.75, adj = 0, line = ifelse(is.na(supp.acmes), 9, 8))
plot(x = NULL, y = NULL, xlim = c(-0.2, 0.2), ylim = c(1,length(supp.nc.diffs)),
     pch = 19, yaxt = "n", bty = "n", xlab = "Eliminated Effect, Support",
     ylab = "", main = "Candidate Party = Republican")
abline(h = 1:length(supp.nc.diffs), col = "grey90")
abline(v = 0, col = "grey60")
points(x = supp.nc.diffs, y = length(supp.nc.diffs):1, pch = 19)
segments(x0 = supp.nc.diff.ci95[1,], x1 = supp.nc.diff.ci95[2,],
         y0 = length(supp.nc.diffs):1)
segments(x0 = supp.nc.diff.ci90[1,], x1 = supp.nc.diff.ci90[2,],
         y0 = length(supp.nc.diffs):1, lwd = 2)
axis(side = 2, at = 1:length(supp.nc.diffs),
     labels = rep("", times = length(supp.nc.diffs)))
mtext(side = 2, at = length(supp.nc.diffs):1, names(supp.nc.diffs), las = 1,
      cex = 0.75, adj = 0, line = ifelse(is.na(supp.nc.diffs), 9, 8))
plot(x = NULL, y = NULL, xlim = c(-0.2, 0.2), ylim = c(1,length(supp.co.diffs)),
     pch = 19, yaxt = "n", bty = "n", xlab = "Eliminated Effect, Support",
     ylab = "", main = "Candidate Party = Democrat")
abline(h = 1:length(supp.co.diffs), col = "grey90")
abline(v = 0, col = "grey60")
points(x = supp.co.diffs, y = length(supp.co.diffs):1, pch = 19)
segments(x0 = supp.co.diff.ci95[1,], x1 = supp.co.diff.ci95[2,],
         y0 = length(supp.co.diffs):1)
segments(x0 = supp.co.diff.ci90[1,], x1 = supp.co.diff.ci90[2,],
         y0 = length(supp.co.diffs):1, lwd = 2)
axis(side = 2, at = 1:length(supp.co.diffs),
     labels = rep("", times = length(supp.co.diffs)))
mtext(side = 2, at = length(supp.co.diffs):1, names(supp.co.diffs), las = 1,
      cex = 0.75, adj = 0, line = ifelse(is.na(supp.co.diffs), 9, 8))
dev.off()

r.only <- 6:9

cairo_pdf(file="conjoint-plot-race.pdf", height = 2.5, width = 6, pointsize = 9, bg = NA)
par(mar = c(4.5, 6.1, 2.1, 0.1), mfrow = c(1,3))
plot(x = NULL, y = NULL, xlim = c(-0.25, 0.25), ylim = range(r.only),
     pch = 19, yaxt = "n", bty = "n", xlab = "Marginal Effects, Support",
     ylab = "", main = "Overall Marginal Effects")
abline(h = 1:length(supp.acmes), col = "grey90")
abline(v = 0, col = "grey60")
points(x = supp.acmes[r.only], y = rev(r.only), pch = 19)
segments(x0 = supp.acme.ci95[1,r.only], x1 = supp.acme.ci95[2,r.only],
         y0 = rev(r.only))
segments(x0 = supp.acme.ci90[1,r.only], x1 = supp.acme.ci90[2,r.only],
         y0 = rev(r.only), lwd = 2)
axis(side = 2, at = r.only,
     labels = rep("", times = length(r.only)))
mtext(side = 2, at = rev(r.only), names(supp.acmes[r.only]), las = 1,
      cex = 0.75, adj = 0, line = ifelse(is.na(supp.acmes[r.only]), 9, 5))
plot(x = NULL, y = NULL, xlim = c(-0.25, 0.25), ylim = range(r.only),
     pch = 19, yaxt = "n", bty = "n", xlab = "Eliminated Effect, Support",
     ylab = "", main = "Nom. Party = Rep")
abline(h = r.only, col = "grey90")
abline(v = 0, col = "grey60")
points(x = supp.nc.diffs[r.only], y = rev(r.only), pch = 19)
segments(x0 = supp.nc.diff.ci95[1, r.only], x1 = supp.nc.diff.ci95[2, r.only],
         y0 = rev(r.only))
segments(x0 = supp.nc.diff.ci90[1, r.only], x1 = supp.nc.diff.ci90[2, r.only],
         y0 = rev(r.only), lwd = 2)
axis(side = 2, at = r.only,
     labels = rep("", times = length(r.only)))
mtext(side = 2, at = rev(r.only), names(supp.nc.diffs[r.only]), las = 1,
      cex = 0.75, adj = 0, line = ifelse(is.na(supp.nc.diffs[r.only]), 9, 5))
plot(x = NULL, y = NULL, xlim = c(-0.25, 0.25), ylim = range(r.only),
     pch = 19, yaxt = "n", bty = "n", xlab = "Eliminated Effect, Support",
     ylab = "", main = "Nom. Party = Dem")
abline(h = r.only, col = "grey90")
abline(v = 0, col = "grey60")
points(x = supp.co.diffs[r.only], y = rev(r.only), pch = 19)
segments(x0 = supp.co.diff.ci95[1, r.only], x1 = supp.co.diff.ci95[2, r.only],
         y0 = rev(r.only))
segments(x0 = supp.co.diff.ci90[1, r.only], x1 = supp.co.diff.ci90[2, r.only],
         y0 = rev(r.only), lwd = 2)
axis(side = 2, at = r.only,
     labels = rep("", times = length(supp.co.diffs[r.only])))
mtext(side = 2, at = rev(r.only), names(supp.co.diffs[r.only]), las = 1,
      cex = 0.75, adj = 0, line = ifelse(is.na(supp.co.diffs[r.only]), 9, 5))
dev.off()
